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Abstract 

We present a simple technique for the computation of coarse-scale steady states of dy- 
namical systems with time scale separation in the form of a "wrapper" around a fine-scale 
simulator. We discuss how this approach alleviates certain problems encountered by com- 
parable existing approaches, and illustrate its use by computing coarse-scale steady states 
of a lattice Boltzmann fine scale code. Interestingly, in the same context of multiple time 
scale problems, the approach can be slightly modified to provide initial conditions on the 
slow manifold with prescribed coarse-scale observables. The approach is based on appro- 
priately designed short bursts of the fine-scale simulator whose results are used to track 
changes in the coarse variables of interest, a core component of the equation-free framework. 
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1. Introduction 

For many problems in science and engineering, the best available model is given at 
a fine-scale level, while we would like to analyze its behavior at a much coarser-scale 
"system level". To bridge the gap between the scale of the available model and the scale 
of interest, one typically attempts to derive a reduced model in terms of an appropriate 
set of variables (the "coarse observables"). In many cases, the derivation of such a reduced 
model hinges on the existence of a low-dimensional, attracting, invariant slow manifold, 
which can be parametrized in terms of these observables. In fine-scale simulations, all 
initial conditions are quickly attracted towards this slow manifold, on which the reduced 
dynamics subsequently evolve. In other words, the remaining fine-scale model variables 
quickly become functionals of (become "slaved to") the observables, and the fine-scale 
model state can be accurately described in terms of the observables only. 
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Although it should in principle always be poss ible to derive a reduced dy nam ical model 
for a s ystem possessing a slow manifold (see, e.g., ISegel and SlemrodI (119891 ) and iRamshaw 
( Il980l )). one may fail to do so, for instance because the closures required in its construction 
are not accurate or even unavailable. F or this class of pr o blems , Kevrekidis and coworkers 
proposed the equation-free framework ( iKevrekidis et al.l (l2003l )). which allows one to per- 
form coarse-scale computations based on appropriately initialized short fine-scale simula- 
tions. (The term "equation-free" emphasizes that the derivation of reduced equations is cir- 
cumvented.) In this paper, we propose and study two related equation- free algorithms that 
are designed to (1) compute stable or unstable coarse-scale steady state solutions, or to (2) 
systematically initializ e a fine-scale model with prescribed value s of th e coarse observable s 
(see, e.g., prior work inlLarn and Goussia (|l994r) 
Gear and KevrekidisI (l2005h : iGear et al.l ( 120051 ) 

( 2002 ): Lorenzl ( 19861) and some further mathematical analysis in lKreiss and Lorenj (Il994j ): 
Zagaris et al. ( 20041 )). As we will see below, the gain in efficiency obtained with these al- 
gorithms stems from the fact that the bulk of the computations are only performed in the 
space of the observables, which is typically low- dimensional compared to the full fine-scale 
variables space. 

The outline of this paper is as follows: In Section |2l we address the computation of 
coarse-scale steady state solutions. We succinctly illustrate certain problems encountered 
by existing, time-stepper based algorithms, and propose a modification that alleviates these 
problems. Section [3] then discusses a variation of the latter approach that helps initialize 
the fine-scale simulator on the slow manifold. Section H] uses a lattice Boltzmann fine scale 
simulator to demonstrate the application of both algorithms; the error characteristics of 
the algorithm is presented in Section [5] and we conclude with a brief discussion. 



2. Computing coarse-scale steady state solutions 

In this section, we outline two existing methods to compute coarse-scale steady state 
solutions, as well as the modification leading to the proposed, third method, with the 
help of a simple model problem. Specifically, we consider the two-dimensional linear time 
integrator 
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as our fine-scale simulator with x being the "coarse observable" . This arises, for instance, 
as the explicit Euler integrator for the system 
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The slow manifold is then y = tan(a)x; note that this manifold can be parametrized by (is 
the graph of a function over) our chosen observable as long as a ^ ±7r/2. Any fine-scale 
initial condition {x, y) is quickly attracted to this slow manifold along a trajectory that (at 
least far enough away from the slow manifold) approximates a line with slope tan(/3). 

2.1. Description of the methods 

The obvious fine-scale steady state solution of ([1]) is (x*, y*) = (0, 0), and this is a stable 
steady state. In this section, we initially illustrate two existing methods whose purpose 
is to approximate the coarse-scale steady state solution x* = 0. These two methods are 



based directly on the concept of the coarse time-stepper (ITheodoropoulos et al.l (j2000[ l: 
Kevrekidis et al.l (120031 )). They approximate the value of x* by solving 
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where $(x, r) denotes a coarse time step over time r with initial condition x. Each such 
coarse time step (of time r) consists of the following three substeps: 

1. lifting, in which an appropriate fine-scale state is constructed according to the value 
of the observable x. 

2. simulation, in which the fine-scale state is evolved over time r, with r large enough 
to allow y to get slaved to x, but small compared to the coarse time scales. 

3. restriction, in which the value of the observable is extracted from the resulting fine- 
scale state. 

If the value of x does not change substantially during the fast transients towards the 
slow manifold, it may suffice to choose arbitrary (say, consistently the same) values for the 
remaining fine-scale variables (in this case, just y) in the lifting step. This type of lifting 
will from now on be called arbitrary lifting, and the resulting method for approximating 
the coarse-scale steady state solution will be called Method 1. 

If the value of x does change substantially during the fast transients towards the slow 
manifold, a more a ccurate initialization (i .e., an initializatio n close r to the slow manifold) 
should be used. In iGear and KevrekidisI ( 120051 ): iGear et al.l (|2005[ ). it was shown that an 
accurate initialization can be obtained with the so-called constrained runs lifting scheme. 
In its simplest form, this scheme determines, for a value x° of the observable, the value 
of y so that dy{x^, y)/dt = 0. The intuitive reason why this condition (or, more generally, 
a condition demanding that the y-time derivative to be bounded) yields a state {x^,y) 
close to the slow manifold is that time differentiation amplifies fast components more than 
slow components, so that, if the time deriva tives are small, th e fast components in the 
remaining fine-scale variables are small. In iGear et al.l ( l2005l ) it was rigorously shown 
that, under certain conditions, the resulting state is indeed a good approximation to the 
desired point on the slow manifold. In practice, it is often convenien t to approximate the 
deriva t ive dy(x^,y)/dt nu merically, e.g. using forward differences. In lGear and Kevrekidis 
( l2005h : lGear et al.l (l2005l ) it was shown that a functional iteration can then, in many cases, 
be used to find the zero of the resulting forward difference condition. If the step size of the 
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functional iteration and the forward difference formula are both equal to the step size of 
the fine-scale simulator At, the functional iteration takes the following form: 

0. Initialize y as well as possible. Then start the iteration 1-3. 

1. Evolve the fine-scale simulator over one time step of size At, starting from y). 

2. Reset the value of x to its original value 

3. If the difference between the current and the previous value of y is smaller than a 
certain tolerance tol: end the iteration. Else: go to 1. 

The iterative scheme above will further be called the constrained runs functional iteration 
(abbreviated CRFI), and the method for approximating the coarse-scale steady state so- 
lution that consists of solving equation (jl]) with y systematically initialized this way will 
be called Method 2. 

In some cases, Method 1 may produce very inaccurate results when the value of x 
changes substantially during the fast transients toward the slow manifold. METHOD 2, 
on the other hand, may find an accurate approximation to the exact solution x* (more 
precisely, the error can be made arbitrarily small by decreasing the value of tol). This can 
easily be seen from the fact that 

• dx/dt (or its finite difference approximation) is zero in the coarse-scale steady state 
solution (this is just equation (jl])). 

• dy/dt (or its finite difference approximation) is zero at the fixed point of the con- 
strained runs functional iteration (because of the nature of this iteration). 

In essence. Method 2 computes the fine-scale steady state solution as a "splitting scheme" , 
by solving the system dy/dt = within each step of an outer solver for the system 
dx/dt = 0. As a result, the computational complexity of Method 2 may be as large 
as that of directly solving the full fine-scale model (which is exactly what we wanted to 
avoid). Moreover, Method 2 may fail when the constrained runs functional iteration does 
not converge to the correct fine-scale state near the slow manifold (the iteration may actu- 
ally be unstable, or it may converge to a solution that does not correspond to a state close 
to the slow manifold). In some cases, these convergence issues may be overcome by using a 



more advanced computational approach such as the one presented in IVandekerckhove et al. 



( 120091 ). yet then the issue of overall computational complexity also remains. 

To cope with the potential accuracy, convergence or efficiency issues, we now propose 
a third method. Method 3, to compute coarse-scale steady state solutions. Instead of 
solving @, this method solves 

^{x,t + t') -^{x,t) = 0, (5) 

in which $(x,r), $(a;,r + r') denote coarse time-steps over times r, r + r' in which we 
use the arbitrary lifting scheme. As before, the value of r should be large enough to 
allow y to get slaved to x, but small compared to the coarse time scales. The variable 
r' represents the time interval of an additional simulation on (or very close to) the slow 
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manifold. If x denotes the solution of ([5]), we expect r) = r + r') (and not x) 
to be a very good approximation to the exact coarse-scale steady state x*, as the finite 
difference approximation of dx/dt based on two points on the slow manifold is then zero. 
The obvious advantages compared to Method 2 are that this method is conceptually 
simpler, that it does not involve the (potentially unstable) constrained runs functional 
iteration, and that we are no longer solving any systems of the form dy{x^,y)/dt = in 
the space of the remaining (fine scale) variables (which may be an advantage when that 
space is large compared to the space of the observables y). 

To summarize: we have outlined three different methods to find coarse-scale steady 
state solutions. In Method 1, we solve ([1]), in which $(x,r) represents a coarse time 
step based on arbitrary lifting. In Method 2, we solve (jlj), in which $(a;, r) represents 
a coarse time step based on the constrained runs functional iteration. In Method 3, we 
solve ([5]), in which <l>(x, r) represents a coarse time step based on arbitrary lifting. 

2.2. Numerical illustration 

We now illustrate the performance of the three methods described above for the model 
problem ([1]) with a = 7r/6 and /3 = — 7r/3. As the smallest eigenvalue of the Jacobian 
matrix of the time integrator ([1]) is A2 = 0.1, it takes about 16 time steps for an 0{1) 
initial condition to reach the slow manifold y = x/Vs up to machine precision. 

Some examples of a single iteration with Method 1 are given in Figured] (top- left). 
The slow manifold is represented by the thick line; the full model steady state solution 
(x*, y*) = (0, 0) on the manifold is indicated by the filled square. For various initial values 
of X, indicated by the filled circles, we perform a time step with a coarse time-stepper 
that is based on the arbitrary lifting scheme x 1— )■ (x, 1/2) and with r = 50 At (remember 
that At is the time step of the fine scale time-stepper ([T])). The fine-scale trajectories are 
represented by the fine lines; note that the slope of these trajectories is approximately 
tan(— 7r/3) = —y/S. The end points of the trajectories are indicated by the open circles; 
these points clearly lie on the slow manifold. The solution of Method 1, the x-coordinate 
of the small open circle, is also shown; the initial condition of the corresponding simulation 
trajectory is encircled as well. For this trajectory, the x-values of the initial and end points 
are equal, as demanded by dl]). In this case, however, the value of x found is 0.719, which 
is clearly different from the exact value x* = 0. As mentioned above, this is due to the 
fact that the value of x changes substantially during the fast transients towards the slow 
manifold (/5 ^ ±7r/2); clearly, the solution found with Method 1 depends on the value 
of T. 

An illustration of Method 2 is given in Figure [T] (right). Again, the slow manifold 
is represented by the thick line and the full model steady state solution {x*,y*) = (0,0) 
is indicated by the filled square. From various initial values of x, indicated by the filled 
circles, we perform the constrained runs functional iteration with tol = 10~^^, after which 
we perform an additional simulation over the time interval r = 50At (here, we only chose 
such a large value of r to make Figure 1 more clear; in practice there is little reason to use 
such a large value of r). The simulation trajectories and the resetting of the observable 
in the constrained runs functional iteration are represented by the fine lines and arrows. 
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Me-thod 1 



Method 2 




Figure 1: An illustration of the three methods used to compute coarse-scale steady states. Shown are the 
slow manifold, fast transients, and computed coarse-scale steady states. See text for further explanation. 



respectively. The end points of the constrained runs functional iteration and the simulation 
trajectory afterwards are indicated by the open circles. We observe that, in this example, 
the constrained runs functional iteration always brings us very close to the slow manifold. 
The final point reached by Method 2 is also shown; the initial condition of the corre- 
sponding simulation trajectory is encircled. As in the case of Method 1, the initial value 
of X is the same as the end value of x, as demanded by (jl]). In this case, however, the value 
of X found is 3.02-10^^^, which, as expected, approximates the exact value x* = up to the 
tolerance tol. In the lower right of this figure, one can see the details of Method 2: first, 
the constrained runs functional iteration is performed until the difference in successive y 
values is less than tol, leaving us at the (blue) triangle; next, we perform a coarse time 
step of length r = 50 At using the ^/-coordinates of this triangle for the lifting, shown as the 
curvy (red) arrow; finally, because the x-coordinate of the (blue) triangle coincides with 
the head of the curvy (red) arrow (meaning that $(x,r) — x = 0), we determine that we 
are at a coarse-scale steady state. 

Due to the fact that the constrained runs lifting brings us very close to the manifold, 
we may even use a smaller value of r than in Method 1; even if r = At, we obtain x = 
1.05-10^^^. It is also worth mentioning that if we had chosen /3 = 7r/4, the constrained runs 
functional iteration would not have converged, as the iteration is then u nstable. For our 
model problem, this can easily be rationalized using geometric arguments (jVandekerckhove 
fl2008h l 

An illustration of Method 3 is given in Figured] (bottom-left). Again, we used the 
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same line styles and markers as before, but now both the end points of the simulation 
trajectories over time r = 50At and r + r' = lOOAt are indicated by the open circles (as 
before, we only chose such a large value of r' to make Figure 1 more clear; in practice 
there is little reason to use such a large value of r'). The solution of Method 3 is also 
shown; the initial condition and the end point of the corresponding simulation trajectory 
are encircled. For this solution, the initial value of x is not the same as the end value of 
x; yet the x-values of the end points of the simulation trajectories after the time interval 
r and r + r' coincide, as demanded for in (|5]). The value of x found is now 1.89 ■ 10~^^, 
which corresponds to the exact value = up to machine precision. Even if r' = At, we 
obtain x = 3.29 ■ 10"^^. 



3. Initializing on a slow manifold 

Remarkably, the modification we presented as Method 3 above to improve the ap- 
proximation of coarse-scale steady state solutions can form the basis of an algorithm for 
appropriately initializing our fine-scale simulators given a desired value of the observable (s). 
In the context of our simple example, this means finding the point y) on the slow man- 
ifold corresponding to some prescribed value of x we denote We already showed how to 
use the constrained runs functional iteration (used in Method 2) for this purpose. The 
fixed point of the constrained runs functional iteration lies close to (but, if x^ is not a 
coarse-scale steady state solution, not exactly on) the slow manifold. More accurate ini- 
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for a ce rtain value of m g N. The larger the value of m, the more 



d™+V(a;°,l/)/dt™+' 

accurate the procedure will be (iGear et al.l (l2005l )). assuming it converges. For these vari- 
ants of the scheme, however, the constrained runs functional iteration is not guaranteed 
to converge to a solution close to the slow manifold, and the computational complexity 
may be unacceptably large. For these reasons, we now propose an alternative initializa- 
tion method that has many similarities with Method 3, but now, instead of computing 
coarse-scale steady state solutions, we compute points lying on the slow manifold for a 
given value of the observable, x^. 

Instead of demanding that the finite difference approximation of the time derivative of 
the observable is zero after a simulation over time r as in ([5]), we now compute x so that 



$(x, r) 



X 



0. 



(6) 



As in Method 3, ^{x, r) denotes a coarse time-step over time r in which we use the 
simple arbitrary lifting scheme. Again, the value of r should be large enough to allow y to 
get slaved to x, but much smaller than the coarse time scales. If x is the solution of ([6]), 
we expect the fine-scale solution {x,y), obtained after simulation over time r starting from 
the arbitrary lifted state corresponding to x, to be a good approximation of the desired 
point on the slow manifold (the simulation step has brought us close to the slow manifold 
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Table 1: Absolute value of the error in the solution of the constrained runs functional iteration and 
InitMan. 
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at the desired value The obvious advantages of this algorithm (which we will refer to 
as InitMan) compared to the constrained runs functional iteration are that this method is 
conceptually simpler, that it does not have the same potential convergence issues, and that 
we are no longer solving a system of equations of the form dy{x^ ,y)/dt = in the space 
of the "remaining variables" (which may be an advantage if that space is large compared 
to the space of the observables) . We still do, of course, require a solver for equation ([H]), 
as would also be required for equation (jlj) above; this might be something as simple as 
a Newt o n or Broyden method, or as sophisticated as a Newton-Krylov GMRES solver 
flKellevI fll995f )l 



To illustrate the performance of the (variants of the) constrained runs functional itera- 
tion and InitMan, we again consider the model problem ([T]) with a = tt/Q and /3 = — 7r/3. 
Using five different values of m for the constrained runs functional iteration, and also 
using InitMan, we approximate the value of y so that {l,y) lies as close as possible 
to the exact point on the slow manifold, (l,l/\/3)- Table [1] shows, for various values 
of m, the error in the solution found by the constrained runs functional iteration with 
tol = 10^^^, and also the error found by InitMan (solved with a Newton iteration). We 
clearly observe that, as the value of m increases, the error decreases by a factor o f about 
(1 - 0.999)/(l - 0.1) ^ 1.11 - 10-3, as expected by theory fIVandekerckhovel fl2008h l The 



errors in the solution found by both the m = 4 constrained runs functional iteration and 
InitMan with r = 15At are at the level of machine precision. 

4. Application to a lattice Boltzmann model 

In this section, we will apply the algorithms we have presented to a lattice Boltzmann 
model (LBM) of a one-dimensional reaction-diffusion system. In Section 14. H we present 
the LBM. In Section we analytically derive two reduced models in terms of two different 
coarse observables. In Section 14.31 a detailed description of the numerical results is given. 

4.I. The lattice Boltzmann model 

An LBM ( Chopard et al. ( 20021 )) describes the evolution of discrete (particle) distribu- 



tion functions fi{xj,tk), which depend on space Xj, time and velocity fj. For our one- 
dimensional model problem, only three values are considered for the velocity {vi = iAx / At, 
with i G {—1, 0, 1}), and each distribution function /j is discretized in space on the domain 
[0, 1] using a grid spacing Ax = 1/N {N lattice intervals) and in time using a time step 
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At. The LBM evolution law for the distributions fi{xj,tk) in the interior of the domain is 

fi{xj+u tk+i) = fi{xj + iAx, tk + At) 

= fi{xj,tk) - uj {fiixj, tk) - ft\xj, tk)) + ^F{p{xj,tk)), 

with i G {-1,0, 1}. 

Diffusive collisions are modeled by the Bhatnagar-Gross-Krook (BG K) collision term 
—uj(fi (xi, tk)—fi'^{xj, tk)) as a relaxation to the local diffusive equilibrium (IQian and Orszag 
fll995h ) 

f!\xj,tk) = -p{xj,tk). (8) 

The parameter w G (0,2) is called the relaxation coefficient and p is the (particle) density 
field, which is defined as the "zeroth" order velocity moment of fi{xj,tk) 

1 

p{Xj,tk) = ^ fi{Xj,tk) = f-liXj,tk) + fo{Xj,tk) + fliXj,tk). 
i=-l 

It follows directly that the BGK diffusive collisions locally conserve density. 

The last term in eq uation (171) models the re a ctions, which a r e assu med to depend only 
on the density field p ( Qian and Orszag ( 1995 ): Dawson et al. ( 1993[ )). In this paper, we 
will use the specific reaction term 



F{p{xj,tk)) = \p{xj,tk) (1 - p{xj,tk)) 



(9) 



in which the parameter A > determines the strength of the reaction "force". Nonlin- 
ear reaction te r ms of the form dHJ a rise naturally i n the fields of heat and mass transfer 
(iDanilov et aD (Il995[ )) or in ecology ( Holmes et al.l (Il994l )). 

At the boundaries, we impose Dirichlet boundary conditions p(0,tfc) = p{^,tk) = by 
assigning the appropriate values to the distribution functions that stream into the domain 
at Xo = and xn = 1- 

Similar to the density p, we can define the momentum and the energy ^ as (a rescaling 
of) the first and the second (or in short, the higher) order moments of fi 



[X 
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tk) = = -f-liXj,tk) + fl{Xj,tk), 



i=-l 
1 



^{Xj,tk) = ^Y1 



f-l{Xj,tk) + fl{Xj,tk) 



i=-l 

Later on we will also use the variable a, which is defined as 
1 

a{xj,tk) = ^{2i'^-l)fi{xj,tk) = f-i{xj,tk)-fo{xj,tk)+fi{xj,tk) = -p{xj,tk)+4:^{xj,tk). 

i=-l 
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4-2. Analytical coarse- graining 

For the LBM described in the previous section, a Chapman- Enskog multiscale expan- 
sion can be used to der i ye an a ccurat e redu ced model for the long-term behavior of the 
system f Chopard et al. ( 2002 ): Succi f 200ll )). (In practice, the equation-free methods 



should of course only be used when such a derivation is not possible. Here, however, the 
analytical derivation provides insight into the problem, which is particularly helpful in 
understanding the performance of the different methods.) In this section, we will derive 
two different reduced models: one in terms of the density p and another in terms of the 
vari able a. For a det a iled d erivation of the following two equations, we refer to appendix 



B of lVandekerckhovd (120081 ). 

If we use the density p as the observable, the reduced model is the partial differential 
equation 



dt dx"^ \ 3u At J dx 

with Dirichlet boundary conditions p{0,t) = p{l,t) = 0. If we consider the variable a to 
be the observable, we obtain the partial differential equation 

d.(.,t) _ ^aV(M) ^ 1^,3^) ^ ('-^^) + 1f(3.) (11) 



dt dx^ 3 ' ' \ 3lu At J dx^ 3 

with Dirichlet boundary conditions cr(0,t) = cT(l,t) = 0. 

As a by-product of the Chapman-Enskog expansion, we find, after dropping the indices 
j and k, and retaining only terms up to second order, that the relation between a and p is 
given by 

, , pixA) 2(2 — u) d'^p(x,t) ^ , n, 
a{x, t) = + \ J ; ^ Ax" + 0{Ax^) 

3 S^uj^ ox^ 1^2^ 



p{x,t) = 3a{x,t) - ^ , ^ ' Ax^ + OiAx""). 



This shows that in a LBM simulation, the value of p is approximately three times as large 
as the value of a, at least after a short initial transient. From this point of view, the choice 
of a as the observable is as natural as the choice of p. 

The fact that the reduced dynamics can be described in terms of p or a only implies 
that the remaining "fine-scale variables" and ^ quickly become functionals of (slaved to) 
p or o". These functionals, which we will call slaving relations, are 

<^(x, t) = - Ax + 0{Ax') (13) 

3uj ox 



e x, t) = -p X, t) - —^L^Ax^ + 0{Ax^) 14 
3 IScj^ ox^ 

in terms of the density p, or 



u ox 
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Figure 2: The full LBM steady state solution, computed in both the (p, 0, ^)- and the (cr, 0, ^)-coordinate 
systems using solid line and circle markers, respectively. This solution is shown in pj a space (left), space 
(middle), and ^ space (right). Dot markers are also present to demonstrate the validity of (fTO|) and 
(left), (dsi) and ([^ (middle), and (Ull) and ^ (right). 



i{x^ t) = a{x, t) + ^^^^A^' + 0{Ax') (16) 

in terms of the variable a. The slaving relations ( |T3l) -(fT^ or ([T5l)-( IT6l) define a slow manifold 
in the LBM phase space, on which the reduced dynamics takes place. 

The validity of the reduced models (fTO|) - (fTT|) and the slaving relations (fT2|) - (fT6|) is 
illustrated in Figured Here, we chose = 100, A = 25 and u = 1.25, and computed the 
full LBM steady state solution in both the (p, 0, ^ )- and the (cr, 0, ^ )-coordinate systems. 
These solutions are shown using the solid line and the circle markers, respectively. To 
simplify the interpretation, we split the solution into the p-, a-, 0- and .^-components, and 
mapped them onto the spatial domain [0, 1]. From these solutions, it is confirmed that the 
value of a is indeed about one third of the value of p. Then, we also computed the steady 
state solutions of ffTOl) and f|TT]) . and added them in the form of the dot markers to Figure 
[2] (left). At the resolution shown here, these reduced solutions are clearly indistinguishable 
from the p- and cr-components of the full LBM solution. Finally, we added to Figure |2] 
(middle and right), also in the form of the dot markers, the profiles of and ^ according to 
the slaving relations (fT3|) - (fT^ (exactly the same results are obtained when using (fT5l) - (fT6l) ). 
Here, we approximated the spatial derivatives of p numerically using finite differences. In 
this case also, the values of and ^ are indistinguishable from the 0- and ^-components of 
the full LBM solution. 

Any LBM initial condition is quickly attracted towards the slow manifold during the 
transient LBM simulation. For u = 1.25, for instance, the slow manifold is reached, up to 
machine precision, in about 25 LBM time steps (the values of and ^ in the simulation 
converge linearly to the slaving relations with convergence factor |1 — u;|, as can be seen 
from d?])). The value of the observable may however change during the transient phase 
towards the slow manifold: 

• Suppose that we start from the initial condition (p, 0, ^) = (po, 0, 0). In terms of the 
distribution functions, this corresponds to (/_i, /o, /i) = (0, po, 0). During the initial 
transient phase, the values of the distribution functions will quickly be redistributed 
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to /o, /i) ~ (po/3, Po/3, Po/'i), or, in terms of the velocity moments, to (p, 0, ~ 
(po, 0, Po/3) (this follows directly from ©-([H]) and fact that the evolution towards the 
slow manifold is very fast compared to the diffusive or reactive time scales). Hence, 
the density p does not change substantially during the transient phase. 

• Suppose that we start from the initial condition (a, 0, ^) = (do, 0,0). In terms of 
the distribution functions, this corresponds to (/-i,/o,/i) = (0,— o"o,0). During 
the initial transient phase, the values of the distribution functions will quickly be 
redistributed to (/_i,/o,/i) ~ (— (Jo/3, — cro/3, — o"o/3), or, in terms of a and the 
higher order velocity moments, to (a, 0, ,^) ~ (— cro/3,0, — cro/3). Hence, the variable 
a does change substantially during the transient phase from ctq to —cro/3. 

As we demonstrated before (and will see below), the fact that the observable changes 
substantially during the transient phase towards the slow manifold may have important 
consequences for the performance of the different methods. 

4.3. Numerical results 

In all numerical experiments reported below, we use D = 1, Ax = 1/N = 1/100 and 
oj = 1.25. According to (ITOD-(IIID, the value of the LBM time step is then At = 2 ■ 10"^ 
For A, we will choose A = 25 or A = 5. For A = 25, the LBM ([7])-([9]) exhibits a nontrivial 
stable steady sta te solu tion; for A = 10, the nontrivial steady state solution is unstable 



( IVandekerckhovel ( 120081 )). These solutions will henceforth be referred to as "the stable 
steady state solution" and "the unstable steady state solution" , respectively. 

To solve the nonlinear systems (jl]), (jS]) or ([6]), which we will write below in more abstract 
form as g{x) = 0, we use the most basic implementation of Newton's method in which we 
estimate the required Jacobian matrices dg/dx{xi) in a column-by-column fashion. Note 
that these Jacobian matrices are typically not very large, as their size is determined by 
the dimension of the coarse subspace; still, the cost of the computational linear algebra 
should be taken into consideration. Column / of dg/dx{xi) is the directional derivative in 
the direction of the l-th unit vector e/, which can be approximated as 

ox e 

with e an appropriate small parameter. The resulting linear systems dg/dx{xi)Ax = 
—g{xi) are then solved by Gaussian elimination. For problems for which the coarse subspace 
is very large, it may be more appropri ate to solve th e linear systems using, for instance, a 
Jacobian- free Newton-Krylov method f Kelley ( 19951 )). The use of these advanced methods 
will not be considered in this article. 

4-3.1. Results using p as the observable 

In this section, we apply our algorithms and the constrained runs functional iteration 
to the LBM, using the density p as the observable. As p does not change substantially 
during the transient phase towards the slow manifold, we expect to obtain good results 
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Table 2: Two- norm of the error ||e||2 = \\pi — p*\\2 in the steady state solution for Method 1 [i = 1), 
Method 2 {i = 2), and Method 3 (i = 3), usmg A = 25 and X — 5. p* is known analytically for each A. 





Method 1 


Method 2 


Method 3 


A 


= 25 (stable) 


1.01 ■ 10-1 


1.07- 10-13 


2.58 • 10-13 


A 


= 5 (unstable) 


4.79 ■ 10-1 


1.47- 10-12 


1.85 ■ 10-12 



with all methods. Note that the results reported here are in p erfect correspondence with 



the theoretical results obtained in lVandekerckhove et al.l (120081 ). in which it was shown that 



the coarse time-stepper under consideration is actually a time integrator for a (slightly) 
modified reaction-diffusion system. 

Coarse-scale steady states: application of Method 1. We approximate the p-component of 
the stable or unstable steady state solution with Method 1, in which we use the arbitrary 
lifting scheme p (-> (p, 0, 0) and set r = 25At. We choose e = y^, with rj ^ 2.22 ■ 10"!^ 
the machine precision. As the initial condition for Newton's method, we take p° = sin(7rx) 
when A = 25 and = — sin(7rx) when A = 5. For both values of A, the nonlinear residual 
reaches the level of machine precision after 5 Newton steps; the resulting solution is called 
Pi. Table |2] shows the two-norm of the error ||e||2 = ||pi — p*\\2 for A = 25 and A = 5 (as 
before, we use p* to denote the p-component of the full LBM steady state solution). The 
error is not large, but it is also not small (the relative error is about 2% when A = 25 and 
about 6% when A = 5). This can be explained by the fact that p does change, although 
very slightly, during the fast transient phase. Note that changing r will also influence the 
approximation of the steady-state solution. 

Coarse-scale steady states: application of METHOD 2. We now approximate the p-component 
of the stable or unstable steady state solution with Method 2, in which we set t = At 
and use the constrained runs lifting scheme with m = and tol = lO-i'*, starting from 
the initial condition (0,0 = (0,0). Again, we choose e = ^^/ rj, and use p° = sin f vrx) o r 



p° = — sin(7rx) as the initial condition for Newton's method. In I Van Leemput et al.l (120081 ) 



it was shown that for m = and p as the observable, the eigenvalues of the constrained 
runs iteration matrix lie on a circle with center point and radius |1 — As a conse- 
quence, the iteration is stable for all values of w G (0, 2) and for u = 1.25 the convergence 
factor is 0.25 (in other words, about 25 iterations are needed to reach the tolerance tol). 
As in Method 1, the nonlinear residual reaches the level of machine precision after 5 
Newton steps for both values of A; the resulting solution is now called p2- Table [2] shows 
the two-norm of the error | |e| I2 = | |p2 ~ P* 1 12 for A = 25 and A = 5. The error is now clearly 
very small. 

Coarse-scale steady states: application 0/ Method 3. We approximate the p-component 
of the stable or unstable steady state solution with Method 3, in which we use the 
arbitrary lifting scheme p (p, 0, 0) and set r = 25At and r' = At. To avoid numeri- 
cal complications due to (nearly) singular Jacobian matrices, we explicitly eliminate the 
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Table 3: Two-norm of the error in the solution of the constrained runs functional iteration and InitMan, 
when p is the observable. 





A = 25 


A = 5 




CRFI, m = 


9.97-10-5 


7.00 ■ 10" 


-5 


CRFI, m = 1 


2.47- 10-^ 


2.59-10" 


-7 


CRFI, m > 2 


oo 


oo 




InitMan 


1.03 ■ 10-15 


1.19 ■ 10" 


15 



boundary conditions p(0) = p(l) = (so the dimension of the Jacobian matrix dg/dx{xi) 
is (A^ — l)x(A^ — 1) instead of (A^ + l)x(A^ + l)). We also set e = ^ (see also the discussion 
in Section [5]). As the initial condition for Newton's method, we again use p° = sin(7rx) or 
p = — sin(7ra;). For both values of A, the nonlinear residual reaches the level of machine 
precision after about 6 or 7 Newton steps; the resulting solution is now called P3. Table 
[2] shows the two-norm of the error ||e||2 = ||p3 — p*||2 = ||'^(P3)'i") ~ P*\\2 for A = 25 and 
A = 5. Again, the error is very small. 

Initialization on the slow manifold: the constrained runs functional iteration ano? InitM AN. 
To test how the constrained runs functional iteration and InitMan perform when we 
attempt to initialize on the slow manifold, we set up the following experiment. For both 
A = 25 and A = 5, we first perform a LBM simulation of 50 steps, starting from the 
(arbitrary) initial condition (p°,0'',^°) = (a;(l — x),x,sm{'Kx)). This provides us with a 
LBM state (p*, 0*, ^*) "on" the slow manifold which is not a coarse scale steady state. Then, 
we use the constrained runs functional iteration or InitMan to approximate the values 
of 0* and ^* corresponding to p*. For the constrained runs functional iteration, we use 
various values of m, set tol = 10-^^ and start from the initial condition (0, ^ = (0, 0). For 
InitMan, we set r = 25 At, start from the initial condition = p*, and again explicitly 
eliminate the boundary conditions p(0) = p(l) = to avoid numerical complications due 
to (nearly) singular Jacobian matrices. We also set e = ^/r/ (see also the discussion in 
Section [5]). 

The results are summarized in Table [31 For both the constrained runs functional 
iteration and InitMan , and for A = 25 and A = 5, we tabulate the two-norm of the 
error ||e||2 = ||(0*,C*) ~ (0*)OI|2 (we use (0*,^*) to denote the solution found by the 
constrained runs functional iteration or InitMan). The solution of the constrained runs 
functional iteration with m = 1 is more accurate than the solution obtained when m = 0, 
but for values of m > 2 the iteration is unstable (some of the eigenvalues of the constrained 
runs iteration matrix are larger than 1 in magnitude). For InitMan, the nonlinear residual 
reaches the level of machine precision after about 3 or 4 Newton steps; the resulting solution 
is called P4. In this case, 0* and ^* are the values of the higher order moments and ^ 
obtained after a LBM simulation over time r starting from (p4,0,0). The solution found 
by InitMan is clearly very accurate. 
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Table 4: Two-norm of the error ||e||2 — ||fTi — cr*||2 in the steady state solution for Method \(i ~ 1), 
Method 2 (i = 2) and Method 3(« — 3), and usmg A = 25 and A = 5. a* is known analytically for each 
A. 





Method 1 


Method 2 


Method 3 


A 


= 25 (stable) 


1.70 


oo 


3.77-10-13 


A 


= 5 (unstable) 


2.69 


oo 


6.16 - 10-13 



4-3.2. Results using a as the observable 

In this section, we apply our algorithms as well as the constrained runs functional 
iteration to the LBM, this time using the variable a as the observable. As a does change 
substantially during the transient phase towards the slow manifold, we expect to obtain 
poor results with METHOD 1, METHOD 2, and the constrained runs functional iteration, 
but good results with Method 3 and InitMan. Tables S] and E] are the analogues of 
Tables [2] and |3l they tabulate the results of the numerical experiments described below. 
Note that the resu lts reported here are i n per fect correspondence with the theoretical 



results obtained in IVandekerckhove et al.l (|2008[ ). in which it was shown that the coarse 



time-stepper under consideration is actually a time integrator for a (slightly) modified 
reaction-diffusion system. 

Coarse-scale steady states: application of Method 1 . We approximate the a-component 
of the stable or unstable steady state solution with Method 1, in which we use the 
arbitrary lifting scheme a i— )■ (a, 0, 0) and set r = 25At. As before, we choose e = y/r]. 
As the initial condition for Newton's method, we use = sin(7rx)/3 when A = 25 or 
(T° = — sin(7rx)/3 when A = 5 (this choice is motivated by the fact that on the slow 
manifold, the value of a is about one third of the value of p; cf. (fT2|) ). For both values of 
A, the nonlinear residual reaches the level of machine precision after 3 Newton steps; the 
resulting solution is called cxi. Table H] shows the two-norm of the error ||e||2 = — cr*\\2 
for A = 25 and A = 5 (as before, we use a* to denote the cx-component of the full LBM 
steady state solution). The error is clearly unacceptable (the iteration converges to a = 0). 

Coarse-scale steady states: application o/ Method 2. We approximate the a-component 
of the stable or unstable steady state solution with Method 2, in which we set t = At 
and use the constrained runs lifting scheme with m = and tol = lO-i'*, starting from 
the initial condition (0,0 = (0,0). Again, we choose e = y/rj, and use (T° = sin(7rx)/3 
or 0"° = — sin(7rx)/3 as the initial condition for Newton's method. For m = 0, u = 1.25 
and A = (we use A = rather than A = 25 or A = 5 as the iteration is then linear), 
and using a as the observable, the eigenvalues of the constrained runs iteration matrix 
lie in (—1.416, —0.25) U (0.25, 1.416). Also for nonzero values of A, the eigenvalues of the 
(varying) iteration matrix fall outside the unit circle. As a consequence, the iteration is 
unstable and METHOD 2 cannot be used; it cannot even be started (see Table H]). 
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Coarse-scale steady states: application 0/ Method 3. We approximate the cr-component 
of the stable or unstable steady state solution with Method 3, in which we use the 
arbitrary lifting scheme a (ct, 0, 0) and set r = 25At and r' = At. To avoid numerical 
complications due to (nearly) singular Jacobian matrices, we again explicitly eliminate 
the boundary conditions cr(0) = cr(l) = and set e = ^ (see also the discussion in 
Section |5]). As the initial condition for Newton's method, we now use cr^ = — sin(7ra;) or 
0"° = sin(7ra;) (this choice is motivated by the fact that, after the initial transient towards 
the slow manifold, the value of a is about minus one third of the value of a", so that we 
then end up near a = sin(7rx)/3 and a = — sin(7rx)/3; cf. the last paragraph in Section 
I4.2p . For both values of A, the nonlinear residual reaches the level of machine precision 
after about 8 (A = 25) or 20 (A = 5) Newton iteration steps; the resulting solution is now 
called as- Table H] shows the two-norm of the error ||e||2 = ||cr3 — o"*||2 = ||$(a3,r) — cr*||2 
for A = 25 and A = 5. Again, the error is very small. 

Initialization on the slow manifold: the constrained runs functional iteration ano? InitMan. 
We now turn again to the problem of initializing on the slow manifold. To test the con- 
strained runs functional iteration and InitMan, we set up the following experiment. For 
both A = 25 and A = 5, we first perform a LBM simulation of 50 steps, starting from the 
(arbitrary) initial condition ((T°,0°,^°) = (x(l — x), x, sin(7rx)). This provides us with a 
LBM state ((T*,0*,^*) "on" the slow manifold. Then, we use the constrained runs func- 
tional iteration or InitMan to approximate the values 0* and ^* corresponding to a*. For 
the constrained runs functional iteration, we use various values of m, set tol = 10"^'^ and 
start from the initial condition (0,^) = (0,0). For InitMan, we set r = 25 At, start from 
the initial condition = —3a* (this choice is again motivated by the fact that, after the 
initial transient towards the slow manifold, the value of a is about minus one third of the 
value of (7°, so that we then end up near a = a*; cf. the last paragraph in Section 14. 2p 
and again explicitly eliminate the boundary conditions a{0) = a{l) = and set e = ^ 
to avoid numerical complications due to (nearly) singular Jacobian matrices (see also the 
discussion in Section E]). 

The results are summarized in Table O For both the constrained runs functional 
iteration and InitMan, and for A = 25 and A = 5, we tabulate the two-norm of the error 
||e||2 = 11(0*)^*) ~ (0*!^*)||2 (as before, we use (0*,^*) to denote the solution found 
by the constrained runs functional iteration or InitMan). As already indicated above, 
the constrained runs functional iteration is always unstable. For InitMan, the nonlinear 
residual reaches the level of machine precision after about 5 or 6 Newton steps; the resulting 
solution is called 0^4. In this case, 0* and ^* are the values of the higher order moments 
and ^ obtained after a LBM simulation over time r starting from (0^4, 0,0). As before, 
the solution found by InitMan is extremely accurate (one cannot expect to do better, in 
fact). 

5. Numerics 

In all of the experiments above, we obtained accurate results with Method 3 and 
InitMan. As we explained, these methods do not suffer from inaccuracies introduced 
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Table 5: Two-norm of the error in the solution of the constrained runs functional iteration and InitMan, 
when cr is the observable. 





A 


= 25 


A = 5 


CRFI, m = 




oo 


oo 


CRFI, m = 1 




oo 


oo 


CRFI, m>2 




oo 


oo 


InitMan 


4.51 


■ 10-14 


1.12 • 10-13 



by arbitrary lifting or from instabilities which sometimes accompany the constrained runs 
functional iteration. In some cases, however, numerical difficulties may also be encoun- 
tered when applying these methods. Let us illustrate this for Method 3 using p as the 
observable. As before, we choose D = 1, Ax = 1/N = 1/100, u = 1.25 and A = 25, use the 
arbitrary lifting scheme p i— )■ [p, 0, 0), set r' = At, use p^ = sin(7rx) as the initial condition 
for Newton's method, and explicitly eliminate the boundary conditions p(0) = p(l) = 0. 
The value of e, however, will now be set equal to e = ^ (with rj the value of machine 
precision) instead of to £ = and we will vary r from to 25At. 

The results are summarized in Table El As expected, the error (compared to the exact 
solution p*) decreases as the value of r/At increases because larger r values bring the fine- 
scale simulator closer to the slow manifold. However, except for very small values of r/At, 
the condition numbers of the Jacobian matrices encountered within Newton's method, 

= \\dg / dx{xi)\ \ ■ \ \{dg/dx{xi))~^\\, also tend to increase. This can be understood by 
realizing that some components of the observable p are also decaying (at a fast pace, albeit 
slower than that of the remaining fine-scale variables), so that more and more relative 
indeterminacy is introduced as the value of r/At increases. In other words, if we perturb 
the exact solution in the direction of a relatively quickly decaying coarse observable, the 
nonlinear residual will remain small as this observable will largely be damped out by the 
time we reach the slow manifold. Since the norm of the Jacobian matrix itself remains 
nearly constant, it is the norm of the inverse of the Jacobian matrix that increases along 
with the condition number. 

As soon as the condition number and norm of the inverse of the Jacobian reach values 
larger than 10^, we observe that Newton's method no longer converges. This can be 
explained as follows. If $(x,r) and $(x,r + r') are both of 0{1), the absolute error 
in g{x) above (remember, we are solving g{x) = 0) is of 0{rj). Due to round-off and 
truncation error in the finite difference approximation of the Jacobian, the absolute error 
in the elements of the Jacobian matrices dg/dx{xi) is then of 0{ri/e + e) = 0{1Q~^). To 
see the influence of this Jacobian matrix perturbation, we can write 

{dg/dx{xi) + 5g){Ax + 6x) = -g{xi), (18) 

in which 5g represents the perturbation matrix (with elements of O{10~^)) and Sx repre- 
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Tabic 6: As a function of r/At, this tabic gives the number of Newton iteration steps required to reach the 
tolerance tol = 10~^^, the two-norm of the nonlinear residual ||$(p3,T + r') — $(p3,r)||2, the two-norm 
of the error ||$(p3,t) — p*\\2, the maximal value of the condition nimiber Ki, and the maximal norm of 
the inverse of the Jacobian \ \{dg/dx{xi))~^\ \ encountered during the Newton iteration. 



r/At 


# iters. 


residual 


error 




maxj \ \{dg / dx{xi 


1 


4 


1.44e-015 


1.73e+000 


2.62e+003 


3.77e+003 


2 


6 


7.05e-016 


5.89C-001 


5.25C+002 


3.77C+003 


3 


5 


7.16e-016 


1.94e-001 


7.34C+002 


3.24C+003 


4 


5 


1.03e-015 


5.95e-002 


2.37e+002 


3.39e+003 


5 


5 


1.59e-015 


1.77e-002 


1.98e+002 


3.38e+003 


6 


5 


9.36e-016 


5.12e-003 


1.70e+002 


3.36e+003 


7 


5 


1.41e-015 


1.45e-003 


4.73e+002 


1.07e+004 


8 


5 


1.30C-015 


4.07C-004 


1.42e+003 


3.59e+004 


9 


5 


1.37e-015 


1.12e-004 


1.53e+002 


4.26e+003 


10 


5 


5.98e-015 


3.08e-005 


8.18e+003 


2.51e+005 


11 


6 


6.64e-016 


8.37e-006 


1.29e+004 


4.29e+005 


12 


5 


5.25C-015 


2.26e-006 


1.88e+003 


6.78e+004 


13 


6 


2.75C-015 


6.06e-007 


2.10e+005 


8.15e+006 


14 


6 


1.18e-015 


1.62e-007 


1.33e-f005 


5.53e+006 


15 


6 


1.16e-015 


4.30e-008 


2.47e+004 


1.09e+006 


16 


8 


2.19e-015 


1.14e-008 


5.15e+005 


2.41e+007 


17 


10 


4.24e-015 


3.01e-009 


1.85e+006 


9.21e+007 


18 


10 


2.10e-015 


7.89e-010 


3.82e+005 


2.00C+007 


19 


9 


NaN 


NaN 


1.51e+008 


4.12C+009 


20 


3 


NaN 


NaN 


6.12e+008 


3.54e+010 


21 


6 


NaN 


NaN 


1.19e+008 


7.20e+009 


22 


6 


NaN 


NaN 


1.22e+009 


7.69e+010 


23 


5 


NaN 


NaN 


8.19C+009 


3.80C+010 


24 


3 


NaN 


NaN 


2.11e+008 


1.45e+010 


25 


3 


NaN 


NaN 


2.00e+009 


1.43e+011 
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sents the resulting perturbation on Ax. If we neglect the SgSx term, we find that 

6x ^ —{dg/dx{xi))^^ ■ 6g ■ Ax. (19) 
It is well known that our (inexact) Newton method converges if 

\\dg/dx{xi){Ax + 6x) + g{xi 



< 1, (20) 



at least if we start sufficiently close to the solution ( Dembo et al. ( 1982I )). Using (|T9l) . this 
becomes 



\\dg/dx{xi){Ax — {dg/dx{xi)) ^ ■ 5g ■ Ax) + g{x-i}\\ \\6g ■ Ax\ 



\9{Xi)\\ \\g{Xi 



< 1- (21) 



If II ■ II denotes the Euclidean vector norm or its induced matrix norm (i.e., the spectral 
norm), it holds that 

\\Sg-^x\\ ^ ll'^^^ll-IIA^^II . . 

= C-D-\\Sg\\-\\idg/dxix,))-'\\ 



\\dg/dx{x,)\\' 



with C,D E [0,1]. (Here, we also used the fact that Aa; = —{dg/dx{xi))^^ ■ g{xi) =^ 
\\Ax\\ = D ■ \ \{dg / dx{xi))^^\ \ ■ ||5f(xj)||.) Since the values C and D are in practice often 
approximately equal to 1, it follows that Newton's method converges if 

INI-ll(a«))-'ll-pJ||L_<, (23) 

at least if we start sufficiently close to the solution. As \\Sg\\ = O{10^^), this im- 
plies that Newton's method is expected to converge if \\{dg / dx{xi))^^\\ < (9(10*). If 
\\{dg / dx{xi))^^\\ > (9(10®), the method is expected to diverge. These theoretical results 
are clearly confirmed in Table El 

Note that for our LBM model problem, the error in the solution can be made as small 
as (9(10~^°) by using the largest value of r/At for which Newton's method converges. If 
the desired level of accuracy cannot be reached due to numerical difficulties, one may try 
the following "trick": increase the value of e, as we did in Section I4l3] when we set e = ^ 
instead of e = y/rj, in an attempt to reduce the error of the finite difference approximation. 
Remember that the finite difference error consists of round-off error and truncation error 
so that the total finite difference error is oi 0{ri/e + e). If the problem is only mildly 
nonlinear, as in the case of our LBM due to the fact that At is small, the constant in the 
0{6) term is small and the finite difference error can be decreased by choosing e > ^Jr]. 
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6. Conclusions 



In this paper, wc have introduced an approach for the computation of coarse-scale 
steady state solutions as well as an approach for initialization on a slow manifold. These 
methods were compared favorably to previously suggested ones: Method 3 and InitMan 
are quick, accurate, and robust, and they bear a striking similarity to each other. We 
demonstrated the use of each of these methods on a lattice Boltzmann model for a reaction- 
diffusion system, compared them to previously suggested methods, and verified our error 
predictions with both numerical results and numerical analysis. These new procedures 
circumvent the need for long, fine-scale simulations to find coarse-scale steady states, or to 
appropriately initialize the fine-scale simulator. 

Our implementation of the numerical methods in this report has been simple and direct, 
in order to clearly illustrate the methods and analyze their sources of error. Indeed, in real- 
world applications, components like the use of higher order derivatives, the reusing of data 
(for example, for two nearby sets of observables, the corresponding fine-scale initializations 
are probably similar), or more intelligent Newton steps may clearly be used to improve 
performance. 
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